stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -3993.9
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 18 -38.4782 0.00113674 0.00127173 1 1 47
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -38.48
## b0 8.80
## b1 1.66
## s 4.15
## m0[1] 18.32
## m0[2] 18.32
## m0[3] 31.00
## m0[4] 46.87
## m0[5] 31.73
## m0[6] 40.71
## m0[7] 20.52
## m0[8] 7.02
## m0[9] 38.01
## m0[10] 37.52
## m0[11] 38.19
## m0[12] 33.96
## m0[13] 38.43
## m0[14] 22.24
## m0[15] 34.45
## m0[16] 17.56
## m0[17] 39.50
## m0[18] 7.40
## m0[19] 29.33
## m0[20] 21.23
## y0[1] 19.40
## y0[2] 16.97
## y0[3] 23.26
## y0[4] 43.19
## y0[5] 31.43
## y0[6] 38.78
## y0[7] 18.25
## y0[8] 11.78
## y0[9] 41.81
## y0[10] 41.95
## y0[11] 33.05
## y0[12] 36.17
## y0[13] 36.94
## y0[14] 19.02
## y0[15] 35.66
## y0[16] 16.87
## y0[17] 33.72
## y0[18] 5.18
## y0[19] 22.95
## y0[20] 22.64
## m1[1] 7.02
## m1[2] 11.45
## m1[3] 15.88
## m1[4] 20.31
## m1[5] 24.73
## m1[6] 29.16
## m1[7] 33.59
## m1[8] 38.02
## m1[9] 42.44
## m1[10] 46.87
## y1[1] 6.69
## y1[2] 11.31
## y1[3] 21.44
## y1[4] 23.48
## y1[5] 22.31
## y1[6] 28.26
## y1[7] 32.13
## y1[8] 44.58
## y1[9] 41.09
## y1[10] 50.81
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -38.76 -38.43 1.35 1.15 -41.46 -37.25 1.00 482 871
## b0 8.73 8.74 2.30 2.35 4.72 12.40 1.02 248 306
## b1 1.67 1.66 0.17 0.16 1.39 1.95 1.02 303 443
## s 4.77 4.64 0.87 0.79 3.58 6.31 1.00 1275 1063
## m0[1] 18.26 18.25 1.53 1.53 15.66 20.66 1.01 276 398
## m0[2] 18.26 18.25 1.53 1.53 15.66 20.66 1.01 276 398
## m0[3] 30.96 30.96 1.14 1.09 29.07 32.77 1.00 1681 1317
## m0[4] 46.85 46.94 2.15 2.06 43.18 50.32 1.01 881 925
## m0[5] 31.69 31.70 1.16 1.12 29.76 33.54 1.00 1935 1315
## m0[6] 40.68 40.72 1.65 1.58 37.89 43.28 1.00 1514 1213
## m0[7] 20.46 20.46 1.39 1.40 18.20 22.65 1.01 301 539
## m0[8] 6.96 6.96 2.46 2.51 2.64 10.88 1.02 248 284
## m0[9] 37.98 38.00 1.46 1.39 35.53 40.27 1.00 1991 1448
## m0[10] 37.49 37.51 1.43 1.36 35.10 39.74 1.00 2060 1517
## m0[11] 38.16 38.18 1.47 1.40 35.69 40.47 1.00 1963 1384
## m0[12] 33.92 33.95 1.24 1.18 31.89 35.85 1.00 2330 1362
## m0[13] 38.40 38.42 1.49 1.42 35.89 40.72 1.00 1928 1384
## m0[14] 22.19 22.18 1.29 1.30 20.06 24.21 1.01 335 668
## m0[15] 34.42 34.44 1.26 1.20 32.35 36.38 1.00 2351 1343
## m0[16] 17.50 17.48 1.59 1.60 14.83 20.00 1.01 269 392
## m0[17] 39.47 39.50 1.56 1.50 36.86 41.90 1.00 1736 1247
## m0[18] 7.33 7.33 2.42 2.47 3.07 11.19 1.02 248 284
## m0[19] 29.29 29.30 1.12 1.10 27.45 31.06 1.00 1009 1254
## m0[20] 21.18 21.17 1.35 1.36 18.95 23.27 1.01 314 600
## y0[1] 18.28 18.41 5.12 5.03 9.91 26.66 1.00 1752 1891
## y0[2] 18.21 18.12 5.18 4.92 9.81 26.80 1.00 1460 1587
## y0[3] 31.01 30.81 4.93 4.67 23.08 38.91 1.00 1662 2013
## y0[4] 46.94 46.80 5.42 5.13 38.01 55.79 1.00 1913 1445
## y0[5] 31.66 31.68 4.94 4.78 23.19 39.70 1.00 2037 1892
## y0[6] 40.89 41.01 5.04 5.00 32.48 49.15 1.00 1822 1636
## y0[7] 20.49 20.64 5.09 4.77 11.87 28.60 1.00 1675 1574
## y0[8] 6.75 6.68 5.41 5.36 -2.08 15.80 1.00 902 1620
## y0[9] 37.97 38.05 5.11 5.01 29.77 45.99 1.00 1864 1972
## y0[10] 37.51 37.55 4.95 4.61 29.43 45.49 1.00 1945 1695
## y0[11] 38.34 38.26 5.08 4.75 30.37 46.57 1.00 2066 1918
## y0[12] 33.95 33.96 5.01 4.80 25.82 42.55 1.00 2124 1889
## y0[13] 38.48 38.36 5.16 4.85 30.11 47.08 1.00 1968 1788
## y0[14] 22.20 22.22 4.99 4.60 13.93 30.51 1.00 1823 1725
## y0[15] 34.23 34.24 5.01 4.86 26.02 42.33 1.00 1990 1947
## y0[16] 17.82 17.75 5.07 4.88 9.96 26.42 1.00 1589 1666
## y0[17] 39.45 39.48 5.06 4.79 31.26 47.62 1.00 2089 1780
## y0[18] 7.03 7.15 5.35 5.08 -1.65 15.63 1.01 689 1354
## y0[19] 29.17 29.31 5.06 4.73 20.52 37.06 1.00 1941 1730
## y0[20] 21.24 21.22 5.07 4.66 13.03 29.38 1.00 1814 1893
## m1[1] 6.96 6.96 2.46 2.51 2.64 10.88 1.02 248 284
## m1[2] 11.39 11.37 2.07 2.14 7.77 14.65 1.02 250 309
## m1[3] 15.82 15.81 1.71 1.73 12.92 18.52 1.02 260 358
## m1[4] 20.25 20.25 1.40 1.40 17.97 22.44 1.01 299 507
## m1[5] 24.69 24.69 1.19 1.18 22.75 26.54 1.01 428 905
## m1[6] 29.12 29.13 1.12 1.10 27.29 30.91 1.00 967 1254
## m1[7] 33.55 33.58 1.22 1.16 31.56 35.48 1.00 2292 1362
## m1[8] 37.98 38.01 1.46 1.39 35.54 40.28 1.00 1990 1448
## m1[9] 42.42 42.47 1.78 1.70 39.36 45.29 1.00 1274 1149
## m1[10] 46.85 46.94 2.15 2.06 43.18 50.32 1.01 881 925
## y1[1] 6.99 6.94 5.40 5.32 -1.91 15.84 1.00 1071 1128
## y1[2] 11.22 11.28 5.37 5.05 2.46 19.87 1.00 1063 1647
## y1[3] 15.87 15.99 5.07 5.07 7.55 24.14 1.00 1427 1553
## y1[4] 20.38 20.35 4.97 4.82 12.23 28.50 1.00 1529 1869
## y1[5] 24.42 24.38 5.00 4.80 16.40 32.71 1.00 1634 1865
## y1[6] 29.09 29.06 4.95 4.68 21.19 37.12 1.00 1910 2057
## y1[7] 33.65 33.72 5.07 4.89 25.47 41.85 1.00 2055 1999
## y1[8] 37.92 38.04 5.09 4.95 29.57 46.02 1.00 2092 1931
## y1[9] 42.41 42.33 5.08 5.07 34.06 50.75 1.00 1869 1790
## y1[10] 46.84 46.85 5.34 5.38 38.32 55.71 1.00 1915 1856
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 0.4 seconds.
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 finished in 0.5 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 finished in 0.6 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.7 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -47.63 -49.21 10.77 9.41 -63.07 -26.53 1.07 54 101
## b0 8.98 8.86 2.37 2.36 5.43 12.85 1.03 143 1209
## b1 1.65 1.65 0.18 0.18 1.37 1.93 1.04 80 1155
## s 3.15 3.24 1.54 1.67 0.70 5.62 1.09 32 26
## sx 1.98 2.00 0.96 0.97 0.48 3.68 1.07 59 70
## x0[1] 3.66 3.81 1.97 2.06 -0.18 6.37 1.05 73 59
## x0[2] 7.02 7.03 1.47 1.62 4.73 9.37 1.03 130 2385
## x0[3] 13.80 13.75 1.23 1.03 11.88 15.88 1.02 1246 1828
## x0[4] 21.65 21.86 1.61 1.61 19.08 24.01 1.04 69 1274
## x0[5] 15.01 14.97 1.49 1.53 12.76 17.47 1.02 192 1391
## x0[6] 20.72 20.54 1.70 1.66 18.29 23.70 1.04 94 417
## x0[7] 5.70 5.79 1.63 1.60 2.87 8.09 1.04 97 448
## x0[8] 0.03 -0.13 1.57 1.58 -2.35 2.68 1.04 77 380
## x0[9] 16.83 16.88 1.31 1.28 14.79 18.91 1.02 168 1896
## x0[10] 16.60 16.67 1.30 1.28 14.52 18.69 1.03 166 1831
## x0[11] 18.39 18.22 1.43 1.27 16.30 20.82 1.02 244 1596
## x0[12] 15.74 15.66 1.26 1.13 13.80 17.87 1.01 862 1525
## x0[13] 17.24 17.33 1.29 1.21 15.24 19.28 1.02 315 2283
## x0[14] 5.02 5.14 2.47 2.85 0.50 8.46 1.06 55 55
## x0[15] 16.42 16.29 1.40 1.31 14.25 18.74 1.02 299 1537
## x0[16] 4.96 5.09 1.38 1.18 2.64 7.05 1.01 597 1326
## x0[17] 20.15 19.93 1.77 1.78 17.72 23.32 1.04 108 788
## x0[18] 0.41 0.26 1.64 1.67 -2.06 3.10 1.04 65 659
## x0[19] 11.12 11.16 1.52 1.62 8.65 13.38 1.03 130 663
## x0[20] 7.06 7.13 1.31 1.10 4.84 9.12 1.01 675 1888
## m0[1] 15.14 14.96 2.97 3.50 10.89 20.09 1.05 61 771
## m0[2] 20.65 20.87 2.47 2.61 16.49 24.25 1.02 230 2675
## m0[3] 31.77 31.90 2.03 1.78 28.42 35.03 1.01 1354 2562
## m0[4] 44.66 44.30 2.56 2.42 40.92 49.20 1.02 336 2502
## m0[5] 33.75 33.88 2.44 2.66 29.73 37.30 1.02 142 2299
## m0[6] 43.13 43.35 2.71 3.00 38.61 46.90 1.03 97 2052
## m0[7] 18.48 18.27 2.53 2.62 14.72 22.73 1.03 93 2255
## m0[8] 9.16 9.48 2.63 2.53 4.56 12.89 1.01 426 2137
## m0[9] 36.75 36.56 2.17 1.92 33.36 40.43 1.01 648 2131
## m0[10] 36.37 36.20 2.09 1.87 33.15 39.84 1.01 600 2420
## m0[11] 39.32 39.52 2.27 2.02 35.48 42.92 1.02 311 2339
## m0[12] 34.96 35.11 2.07 1.91 31.40 38.15 1.02 1116 2676
## m0[13] 37.43 37.30 2.11 1.84 34.09 41.00 1.01 1135 2405
## m0[14] 17.37 17.28 3.83 5.02 12.06 23.42 1.07 46 329
## m0[15] 36.06 36.24 2.27 2.26 32.31 39.45 1.01 331 2927
## m0[16] 17.26 17.27 2.14 1.80 13.81 20.76 1.02 2866 2173
## m0[17] 42.18 42.37 2.77 3.17 37.59 46.10 1.03 90 2191
## m0[18] 9.79 10.13 2.74 2.70 5.11 13.58 1.02 278 1915
## m0[19] 27.37 27.26 2.44 2.57 23.80 31.40 1.03 118 2156
## m0[20] 20.71 20.61 2.08 1.75 17.26 24.16 1.02 1848 2395
## y0[1] 15.18 14.23 4.67 4.23 8.90 23.66 1.02 253 2565
## y0[2] 20.57 21.28 4.30 3.48 12.95 26.79 1.02 1027 2648
## y0[3] 31.83 32.11 3.96 3.23 25.00 38.20 1.02 3404 2528
## y0[4] 44.61 43.98 4.26 3.47 38.35 52.04 1.03 1576 2989
## y0[5] 33.82 34.45 4.27 3.48 26.11 39.94 1.01 951 2539
## y0[6] 43.12 43.84 4.42 3.59 35.17 49.38 1.01 456 2510
## y0[7] 18.53 17.94 4.31 3.47 12.11 26.21 1.02 502 2915
## y0[8] 9.14 9.70 4.34 3.55 1.32 15.72 1.01 1650 2552
## y0[9] 36.69 36.42 4.07 3.19 30.38 43.95 1.02 3185 2333
## y0[10] 36.27 35.95 4.19 3.44 29.87 43.39 1.02 3633 2931
## y0[11] 39.37 39.74 4.21 3.31 31.96 45.87 1.02 2197 2842
## y0[12] 34.97 35.27 4.01 3.15 28.11 41.29 1.02 2681 2709
## y0[13] 37.36 37.08 4.12 3.30 30.95 44.33 1.03 2693 2712
## y0[14] 17.40 16.47 5.14 5.25 10.81 26.89 1.04 84 1993
## y0[15] 35.99 36.50 4.28 3.44 28.40 42.48 1.02 1241 2159
## y0[16] 17.33 17.24 4.13 3.24 10.59 24.29 1.03 3836 2231
## y0[17] 42.11 42.83 4.51 3.91 33.75 48.34 1.01 384 2462
## y0[18] 9.84 10.45 4.51 3.61 1.75 16.44 1.01 1202 3066
## y0[19] 27.46 26.83 4.26 3.56 21.19 35.23 1.01 754 2806
## y0[20] 20.70 20.56 4.11 3.24 14.00 27.50 1.03 3705 3174
## m1[1] 7.22 7.09 2.54 2.53 3.42 11.33 1.03 131 1268
## m1[2] 11.62 11.51 2.11 2.11 8.42 15.06 1.03 166 1368
## m1[3] 16.02 15.93 1.71 1.72 13.42 18.81 1.02 224 1370
## m1[4] 20.42 20.36 1.37 1.32 18.29 22.63 1.01 493 1576
## m1[5] 24.82 24.81 1.12 1.04 22.94 26.56 1.02 1269 2009
## m1[6] 29.22 29.27 1.05 0.95 27.44 30.86 1.01 1824 2137
## m1[7] 33.62 33.61 1.19 1.20 31.67 35.48 1.01 401 2122
## m1[8] 38.02 38.00 1.48 1.59 35.61 40.28 1.03 138 2036
## m1[9] 42.42 42.37 1.85 2.00 39.47 45.26 1.03 106 1685
## m1[10] 46.82 46.76 2.26 2.41 43.31 50.40 1.04 94 959
## x1[1] -1.11 -1.09 2.24 1.67 -4.92 2.69 1.02 4230 2272
## x1[2] 1.55 1.55 2.17 1.61 -2.05 5.17 1.02 3986 1588
## x1[3] 4.24 4.24 2.20 1.61 0.59 7.76 1.02 3708 2226
## x1[4] 7.00 6.99 2.18 1.67 3.50 10.64 1.02 4243 3124
## x1[5] 9.57 9.60 2.15 1.60 6.01 13.08 1.02 4066 2474
## x1[6] 12.31 12.26 2.25 1.67 8.75 16.05 1.02 3840 2435
## x1[7] 14.85 14.87 2.13 1.65 11.30 18.34 1.02 3955 2330
## x1[8] 17.52 17.56 2.24 1.66 13.70 21.16 1.02 4001 2018
## x1[9] 20.20 20.20 2.21 1.68 16.67 23.84 1.02 3772 2267
## x1[10] 22.86 22.87 2.27 1.64 19.18 26.49 1.02 3524 2260
## y1[1] 7.26 7.10 4.35 4.14 0.31 14.28 1.01 749 2586
## y1[2] 11.68 11.54 4.06 3.65 5.13 18.27 1.01 944 2129
## y1[3] 15.95 15.90 3.95 3.29 9.58 22.41 1.01 1321 2520
## y1[4] 20.38 20.31 3.80 2.94 14.31 26.79 1.02 3401 2765
## y1[5] 24.78 24.87 3.69 2.85 18.53 30.70 1.03 3594 2582
## y1[6] 29.24 29.30 3.69 2.83 22.96 35.15 1.03 3755 2313
## y1[7] 33.56 33.75 3.76 2.81 27.12 39.48 1.02 2704 2455
## y1[8] 38.08 38.12 3.87 3.17 31.78 44.49 1.01 1790 3036
## y1[9] 42.34 42.40 3.96 3.44 35.52 48.57 1.01 916 2423
## y1[10] 46.82 46.96 4.20 3.91 40.05 53.55 1.01 515 2761
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -60267.9
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 26 -33.7153 0.000295205 0.000111983 1 1 71
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -33.72
## b0 6.00
## b1 1.86
## s 3.27
## m0[1] 11.57
## m0[2] 21.43
## m0[3] 8.68
## m0[4] 19.77
## m0[5] 13.48
## m0[6] 13.99
## m0[7] 17.35
## m0[8] 13.47
## m0[9] 19.21
## m0[10] 7.15
## m0[11] 9.82
## m0[12] 6.26
## m0[13] 17.35
## m0[14] 15.68
## m0[15] 20.82
## m0[16] 13.73
## m0[17] 20.82
## m0[18] 17.98
## m0[19] 19.16
## m0[20] 20.26
## y0[1] 12.23
## y0[2] 24.51
## y0[3] 10.42
## y0[4] 24.35
## y0[5] 12.97
## y0[6] 18.25
## y0[7] 18.06
## y0[8] 12.42
## y0[9] 20.31
## y0[10] 7.74
## y0[11] 10.88
## y0[12] 1.30
## y0[13] 15.88
## y0[14] 18.45
## y0[15] 24.94
## y0[16] 22.56
## y0[17] 19.99
## y0[18] 19.73
## y0[19] 19.53
## y0[20] 18.98
## m1[1] 6.26
## m1[2] 7.95
## m1[3] 9.63
## m1[4] 11.32
## m1[5] 13.00
## m1[6] 14.69
## m1[7] 16.37
## m1[8] 18.06
## m1[9] 19.74
## m1[10] 21.43
## y1[1] 5.93
## y1[2] 8.75
## y1[3] 6.32
## y1[4] 9.92
## y1[5] 16.23
## y1[6] 13.25
## y1[7] 12.10
## y1[8] 19.00
## y1[9] 15.78
## y1[10] 19.57
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -34.08 -33.75 1.27 1.05 -36.60 -32.69 1.00 650 951
## b0 6.10 6.03 1.76 1.76 3.35 8.96 1.01 394 451
## b1 1.84 1.84 0.32 0.31 1.31 2.36 1.01 474 508
## s 3.71 3.62 0.68 0.61 2.76 4.96 1.00 1268 1136
## m0[1] 11.62 11.62 1.02 0.99 9.96 13.33 1.01 456 692
## m0[2] 21.39 21.38 1.37 1.27 19.11 23.67 1.00 1030 1048
## m0[3] 8.76 8.73 1.37 1.36 6.63 11.00 1.01 400 472
## m0[4] 19.75 19.76 1.15 1.03 17.82 21.68 1.00 1292 1152
## m0[5] 13.51 13.52 0.87 0.83 12.08 14.96 1.01 604 1042
## m0[6] 14.02 14.02 0.84 0.79 12.62 15.41 1.01 678 1157
## m0[7] 17.35 17.37 0.91 0.83 15.83 18.84 1.00 1622 1339
## m0[8] 13.51 13.52 0.87 0.83 12.08 14.96 1.01 603 1042
## m0[9] 19.19 19.20 1.09 0.97 17.38 21.01 1.00 1412 1147
## m0[10] 7.24 7.20 1.59 1.58 4.78 9.82 1.01 394 472
## m0[11] 9.89 9.87 1.22 1.22 7.97 11.87 1.01 412 575
## m0[12] 6.36 6.31 1.72 1.72 3.67 9.15 1.01 393 438
## m0[13] 17.35 17.38 0.91 0.83 15.83 18.85 1.00 1622 1339
## m0[14] 15.70 15.71 0.83 0.77 14.30 17.04 1.00 1099 1304
## m0[15] 20.79 20.78 1.29 1.17 18.65 22.93 1.00 1116 1151
## m0[16] 13.76 13.76 0.85 0.81 12.35 15.19 1.01 638 1090
## m0[17] 20.79 20.79 1.29 1.17 18.65 22.93 1.00 1116 1151
## m0[18] 17.97 18.00 0.96 0.87 16.36 19.56 1.00 1618 1297
## m0[19] 19.14 19.15 1.08 0.97 17.34 20.96 1.00 1422 1147
## m0[20] 20.23 20.24 1.21 1.10 18.19 22.26 1.00 1204 1151
## y0[1] 11.70 11.68 3.95 3.83 5.30 18.30 1.00 1797 1864
## y0[2] 21.42 21.34 3.95 3.70 15.14 28.01 1.00 1739 1733
## y0[3] 8.81 8.75 4.04 3.85 2.13 15.62 1.00 1378 2084
## y0[4] 19.68 19.71 3.92 3.83 13.26 25.85 1.00 1982 1969
## y0[5] 13.52 13.63 3.84 3.80 6.99 19.63 1.00 1859 1717
## y0[6] 14.13 14.07 3.80 3.78 8.18 20.36 1.00 1521 1762
## y0[7] 17.21 17.26 3.78 3.64 10.67 23.42 1.00 2137 1970
## y0[8] 13.39 13.52 3.91 3.79 6.80 19.59 1.00 1805 1670
## y0[9] 19.15 19.11 4.05 3.95 12.66 25.75 1.00 2176 2058
## y0[10] 7.29 7.31 4.07 3.89 0.64 13.89 1.00 1194 1570
## y0[11] 9.89 9.88 4.01 3.73 3.22 16.73 1.00 1330 1843
## y0[12] 6.26 6.21 4.15 4.02 -0.42 12.98 1.00 1220 1416
## y0[13] 17.16 17.11 3.94 3.81 10.64 23.68 1.00 2149 1975
## y0[14] 15.70 15.65 3.86 3.65 9.44 22.19 1.00 2003 1949
## y0[15] 20.86 20.94 4.01 3.91 14.43 27.53 1.00 2020 1917
## y0[16] 13.81 13.73 4.02 3.94 7.04 20.37 1.00 1913 1362
## y0[17] 20.76 20.84 4.04 3.82 14.11 27.43 1.00 2064 1892
## y0[18] 18.02 18.04 3.93 3.89 11.41 24.32 1.00 2005 1770
## y0[19] 18.99 19.01 4.07 3.86 12.25 25.48 1.00 1987 1955
## y0[20] 20.23 20.29 3.87 3.75 13.88 26.53 1.00 1735 1857
## m1[1] 6.36 6.31 1.72 1.72 3.67 9.15 1.01 393 438
## m1[2] 8.03 8.01 1.48 1.48 5.75 10.42 1.01 397 473
## m1[3] 9.70 9.68 1.25 1.26 7.75 11.74 1.01 409 569
## m1[4] 11.37 11.37 1.05 1.02 9.68 13.11 1.01 447 697
## m1[5] 13.04 13.05 0.90 0.87 11.58 14.56 1.01 550 938
## m1[6] 14.71 14.72 0.82 0.78 13.36 16.09 1.00 805 1228
## m1[7] 16.38 16.41 0.85 0.79 14.97 17.77 1.00 1368 1319
## m1[8] 18.05 18.08 0.97 0.88 16.42 19.65 1.00 1605 1297
## m1[9] 19.72 19.73 1.15 1.03 17.80 21.65 1.00 1299 1152
## m1[10] 21.39 21.38 1.37 1.27 19.11 23.67 1.00 1030 1048
## y1[1] 6.40 6.34 4.16 3.96 -0.32 13.27 1.01 1206 1380
## y1[2] 7.97 7.98 4.20 4.08 1.22 14.93 1.00 1369 1959
## y1[3] 9.67 9.66 4.01 3.76 3.11 16.31 1.00 1531 1689
## y1[4] 11.43 11.56 3.94 3.84 4.88 17.95 1.00 1285 1835
## y1[5] 12.95 12.91 3.80 3.60 6.72 19.22 1.00 1803 1812
## y1[6] 14.85 14.81 3.85 3.74 8.56 21.26 1.00 2093 1891
## y1[7] 16.42 16.40 3.74 3.62 10.20 22.37 1.00 1891 1798
## y1[8] 18.01 18.05 3.91 3.75 11.80 24.32 1.00 1861 1922
## y1[9] 19.76 19.72 3.95 3.78 13.29 26.12 1.00 1891 1991
## y1[10] 21.37 21.41 3.93 3.72 14.95 27.85 1.00 1925 1672
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -123.73
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 21 -14.5206 0.00175654 0.000469138 0.9852 0.9852 31
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -14.52
## b0 4.33
## b1 2.04
## s 0.67
## m0[1] 10.44
## m0[2] 21.25
## m0[3] 7.28
## m0[4] 19.43
## m0[5] 12.53
## m0[6] 13.09
## m0[7] 16.77
## m0[8] 12.53
## m0[9] 18.81
## m0[10] 5.59
## m0[11] 8.52
## m0[12] 4.62
## m0[13] 16.78
## m0[14] 14.95
## m0[15] 20.58
## m0[16] 12.81
## m0[17] 20.58
## m0[18] 17.46
## m0[19] 18.76
## m0[20] 19.97
## y0[1] 10.53
## y0[2] 20.86
## y0[3] 8.25
## y0[4] 19.23
## y0[5] 12.43
## y0[6] 13.62
## y0[7] 16.10
## y0[8] 12.95
## y0[9] 19.39
## y0[10] 7.30
## y0[11] 6.65
## y0[12] 4.40
## y0[13] 19.75
## y0[14] 24.10
## y0[15] 20.91
## y0[16] 13.48
## y0[17] 20.55
## y0[18] 17.68
## y0[19] 18.38
## y0[20] 22.28
## m1[1] 4.62
## m1[2] 6.47
## m1[3] 8.32
## m1[4] 10.16
## m1[5] 12.01
## m1[6] 13.86
## m1[7] 15.71
## m1[8] 17.55
## m1[9] 19.40
## m1[10] 21.25
## y1[1] -13.47
## y1[2] 7.75
## y1[3] 4.61
## y1[4] 9.69
## y1[5] 12.94
## y1[6] 15.74
## y1[7] 16.48
## y1[8] 15.49
## y1[9] 19.11
## y1[10] 18.78
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.46 -16.12 1.27 1.03 -19.12 -15.07 1.01 612 686
## b0 3.96 4.04 0.65 0.64 2.77 4.87 1.01 353 433
## b1 2.09 2.08 0.11 0.10 1.92 2.28 1.01 409 513
## s 0.81 0.77 0.24 0.21 0.47 1.23 1.00 752 660
## m0[1] 10.22 10.26 0.39 0.39 9.53 10.79 1.01 388 557
## m0[2] 21.31 21.31 0.41 0.38 20.63 21.97 1.00 1096 1025
## m0[3] 6.98 7.04 0.52 0.51 6.04 7.71 1.01 357 487
## m0[4] 19.44 19.44 0.34 0.32 18.87 20.01 1.00 1484 1155
## m0[5] 12.37 12.38 0.32 0.32 11.82 12.86 1.01 461 657
## m0[6] 12.94 12.95 0.31 0.31 12.42 13.44 1.01 500 664
## m0[7] 16.72 16.71 0.28 0.27 16.26 17.21 1.00 1466 1157
## m0[8] 12.36 12.38 0.32 0.32 11.81 12.86 1.01 460 657
## m0[9] 18.81 18.80 0.33 0.30 18.26 19.35 1.00 1615 1232
## m0[10] 5.25 5.33 0.60 0.58 4.16 6.08 1.01 354 474
## m0[11] 8.25 8.31 0.46 0.46 7.42 8.92 1.01 363 528
## m0[12] 4.26 4.34 0.64 0.62 3.08 5.15 1.01 353 438
## m0[13] 16.72 16.71 0.28 0.27 16.26 17.21 1.00 1467 1157
## m0[14] 14.85 14.84 0.28 0.28 14.39 15.33 1.01 792 1264
## m0[15] 20.62 20.63 0.38 0.35 19.98 21.25 1.00 1225 1005
## m0[16] 12.65 12.66 0.31 0.31 12.11 13.14 1.01 480 698
## m0[17] 20.62 20.63 0.38 0.35 19.98 21.25 1.00 1225 1005
## m0[18] 17.42 17.41 0.29 0.28 16.95 17.93 1.00 1644 1192
## m0[19] 18.75 18.75 0.32 0.30 18.21 19.30 1.00 1625 1232
## m0[20] 19.99 19.99 0.36 0.33 19.39 20.57 1.00 1356 1072
## y0[1] 9.71 10.21 28.05 1.34 3.69 16.25 1.00 1851 1929
## y0[2] 20.36 21.27 49.97 1.31 15.45 25.83 1.00 1686 1587
## y0[3] 6.28 6.97 32.37 1.40 2.15 12.18 1.00 1641 2012
## y0[4] 20.92 19.45 35.11 1.22 14.72 24.54 1.00 1992 1872
## y0[5] 12.23 12.37 19.04 1.21 7.74 17.49 1.00 1633 1707
## y0[6] 12.88 12.92 46.89 1.21 8.06 17.92 1.00 1914 1857
## y0[7] 17.18 16.73 19.70 1.22 12.04 22.59 1.00 2077 1973
## y0[8] 12.07 12.35 18.45 1.20 7.57 16.81 1.00 1950 1929
## y0[9] 19.03 18.83 22.33 1.25 13.97 24.27 1.00 2225 1922
## y0[10] 5.78 5.30 28.75 1.35 -0.74 9.92 1.00 1446 1931
## y0[11] 8.76 8.25 20.27 1.28 3.08 12.48 1.00 1576 1954
## y0[12] 3.47 4.27 22.70 1.45 -0.77 9.75 1.00 1574 1963
## y0[13] 14.83 16.72 102.97 1.20 11.72 20.90 1.00 2058 1960
## y0[14] 22.15 14.80 278.41 1.18 10.24 20.07 1.00 1562 1763
## y0[15] 21.28 20.59 33.42 1.30 15.59 25.97 1.00 1871 2024
## y0[16] 12.31 12.70 31.48 1.19 7.01 17.85 1.00 1933 1828
## y0[17] 15.48 20.61 258.70 1.36 15.04 26.11 1.00 1948 1859
## y0[18] 18.33 17.37 48.11 1.30 11.96 22.64 1.00 2147 1974
## y0[19] 16.65 18.74 108.24 1.20 13.48 23.16 1.00 1967 1861
## y0[20] 20.10 19.96 26.40 1.28 15.30 26.17 1.00 1922 1927
## m1[1] 4.26 4.34 0.64 0.62 3.08 5.15 1.01 353 438
## m1[2] 6.15 6.23 0.55 0.54 5.14 6.92 1.01 355 484
## m1[3] 8.04 8.11 0.47 0.47 7.20 8.72 1.01 362 507
## m1[4] 9.94 9.99 0.40 0.40 9.23 10.52 1.01 382 569
## m1[5] 11.83 11.86 0.34 0.34 11.25 12.35 1.01 436 626
## m1[6] 13.73 13.73 0.29 0.29 13.25 14.21 1.01 583 977
## m1[7] 15.62 15.61 0.28 0.27 15.17 16.11 1.01 1033 1280
## m1[8] 17.52 17.51 0.30 0.28 17.04 18.02 1.00 1659 1122
## m1[9] 19.41 19.41 0.34 0.32 18.84 19.97 1.00 1489 1127
## m1[10] 21.31 21.31 0.41 0.38 20.63 21.97 1.00 1096 1025
## y1[1] 4.88 4.25 15.43 1.48 0.13 9.23 1.00 1223 1960
## y1[2] 3.27 6.13 121.80 1.37 1.28 11.12 1.00 1255 1943
## y1[3] 6.97 8.08 62.10 1.35 2.96 13.20 1.00 1551 2031
## y1[4] 11.30 9.97 40.66 1.32 4.54 15.12 1.00 1539 1915
## y1[5] 12.60 11.84 19.61 1.20 6.92 16.97 1.01 1280 1595
## y1[6] 14.87 13.75 57.63 1.22 8.61 18.83 1.00 1771 1979
## y1[7] 15.24 15.66 32.76 1.28 10.62 21.18 1.00 2010 1788
## y1[8] 18.06 17.48 43.64 1.18 11.90 22.49 1.00 1890 1931
## y1[9] 19.62 19.39 9.69 1.24 15.14 24.86 1.00 2063 1973
## y1[10] 21.49 21.29 22.72 1.24 16.13 26.62 1.00 1828 2004
censored data
objective variable has NA
(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)ex11-0.stan
data{
int N0;
array[N0] vector[2] xy;
int N1;
vector[N1] x1;
}
parameters{
vector[2] m;
cov_matrix[2] s;
}
model{
target+=multi_normal_lpdf(xy | m, s);
x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
vector[2] xy1;
xy1=multi_normal_rng(m,s);
real cr;
cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)## [1] 0.797qplot(x,y)
L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)## [1] 0.456qplot(x0,y0)
mdl=cmdstan_model('./ex11-0.stan')
data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -6228.25
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 56 -105.906 0.000736808 0.00147673 1 1 66
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -105.91
## m[1] 4.26
## m[2] 24.33
## s[1,1] 5.20
## s[2,1] 11.48
## s[1,2] 11.48
## s[2,2] 47.82
## xy1[1] 4.13
## xy1[2] 22.59
## cr 0.73mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.3 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -101.69 -101.31 1.79 1.61 -105.23 -99.50 1.00 582 771
## m[1] 4.25 4.23 0.50 0.48 3.48 5.09 1.00 887 861
## m[2] 24.49 24.69 2.97 2.76 19.34 29.04 1.02 374 445
## s[1,1] 6.80 6.39 2.17 1.84 4.14 11.10 1.00 1468 1697
## s[2,1] 14.41 13.52 9.67 8.97 0.77 31.71 1.01 374 479
## s[1,2] 14.41 13.52 9.67 8.97 0.77 31.71 1.01 374 479
## s[2,2] 76.64 64.30 44.95 34.81 28.01 163.32 1.01 443 525
## xy1[1] 4.26 4.22 2.63 2.52 -0.11 8.64 1.00 2016 1955
## xy1[2] 24.61 24.66 9.17 8.34 9.66 39.38 1.00 1305 1309
## cr 0.61 0.70 0.27 0.20 0.05 0.89 1.01 406 491
xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])## [1] 0.607qplot(xy[,,1],xy[,,2])
objective variable is censored
y i=1-N, y~N(m,s)
actual ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11-1.stan
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionmdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -2060.46
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 22 -11.1021 0.000200091 0.000411365 0.962 0.962 38
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -11.10
## b0 12.87
## b1 2.04
## s 2.08
## m0[1] 23.18
## m0[2] 20.68
## m0[3] 21.08
## m0[4] 28.17
## m0[5] 18.08
## m0[6] 27.39
## m0[7] 18.13
## m0[8] 18.88
## m0[9] 24.07
## y0[1] 23.12
## y0[2] 21.31
## y0[3] 21.43
## y0[4] 27.30
## y0[5] 23.80
## y0[6] 28.29
## y0[7] 19.74
## y0[8] 17.98
## y0[9] 28.01
## m1[1] 12.98
## m1[2] 14.93
## m1[3] 16.88
## m1[4] 18.83
## m1[5] 20.79
## m1[6] 22.74
## m1[7] 24.69
## m1[8] 26.64
## m1[9] 28.59
## m1[10] 30.55
## y1[1] 12.25
## y1[2] 13.97
## y1[3] 16.55
## y1[4] 18.01
## y1[5] 21.85
## y1[6] 19.63
## y1[7] 28.27
## y1[8] 25.39
## y1[9] 30.85
## y1[10] 29.58mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.15 -11.73 1.50 1.25 -15.07 -10.53 1.00 395 441
## b0 13.08 12.98 2.79 2.50 8.65 17.94 1.01 258 295
## b1 2.00 2.02 0.56 0.49 1.07 2.91 1.01 296 380
## s 2.90 2.65 1.05 0.76 1.74 4.91 1.00 700 880
## m0[1] 23.23 23.23 1.03 0.91 21.55 24.94 1.00 1539 1106
## m0[2] 20.78 20.74 1.10 0.98 19.01 22.63 1.00 552 738
## m0[3] 21.17 21.14 1.06 0.95 19.48 22.94 1.00 633 706
## m0[4] 28.15 28.15 1.90 1.75 25.12 31.22 1.00 672 878
## m0[5] 18.22 18.17 1.54 1.32 15.74 20.87 1.01 306 361
## m0[6] 27.38 27.40 1.72 1.58 24.67 30.14 1.00 754 933
## m0[7] 18.26 18.21 1.53 1.32 15.80 20.89 1.01 307 361
## m0[8] 19.00 18.94 1.38 1.17 16.82 21.34 1.00 341 436
## m0[9] 24.11 24.12 1.11 0.98 22.30 25.95 1.00 1544 1093
## y0[1] 23.28 23.31 3.11 2.59 18.27 28.11 1.00 1841 1611
## y0[2] 20.73 20.69 3.31 2.91 15.67 26.05 1.00 1602 1383
## y0[3] 21.15 21.06 3.19 2.67 16.26 26.41 1.00 1821 1531
## y0[4] 28.08 28.06 3.72 3.12 21.82 34.24 1.00 1126 1191
## y0[5] 18.20 18.14 3.50 3.07 12.87 23.90 1.00 1013 1320
## y0[6] 27.33 27.28 3.40 2.95 21.89 32.96 1.00 1460 1488
## y0[7] 18.32 18.21 3.52 3.05 12.97 23.81 1.00 914 1216
## y0[8] 19.02 19.04 3.53 2.95 13.61 24.66 1.00 1028 1373
## y0[9] 24.12 24.11 3.22 2.83 19.05 29.35 1.00 1940 1736
## m1[1] 13.19 13.10 2.77 2.49 8.80 18.02 1.01 258 295
## m1[2] 15.11 15.04 2.27 2.02 11.57 19.07 1.01 264 312
## m1[3] 17.03 16.98 1.81 1.56 14.25 20.13 1.01 280 368
## m1[4] 18.96 18.90 1.39 1.19 16.76 21.32 1.00 338 433
## m1[5] 20.88 20.85 1.09 0.97 19.13 22.68 1.00 571 736
## m1[6] 22.80 22.80 1.01 0.88 21.18 24.48 1.00 1363 1221
## m1[7] 24.72 24.72 1.20 1.04 22.79 26.69 1.00 1355 1173
## m1[8] 26.64 26.65 1.55 1.40 24.16 29.07 1.00 859 992
## m1[9] 28.57 28.57 2.00 1.86 25.34 31.83 1.00 639 838
## m1[10] 30.49 30.50 2.48 2.28 26.42 34.56 1.00 543 760
## y1[1] 13.11 13.13 4.06 3.61 6.44 19.86 1.00 704 956
## y1[2] 15.13 15.08 3.87 3.37 8.96 21.71 1.00 790 1112
## y1[3] 16.90 16.81 3.59 3.16 10.98 22.71 1.00 818 1071
## y1[4] 18.96 19.02 3.39 2.96 13.57 24.40 1.00 1243 1074
## y1[5] 20.88 20.87 3.28 2.79 15.95 26.06 1.00 1776 1736
## y1[6] 22.83 22.84 3.23 2.82 17.90 28.00 1.00 1784 1786
## y1[7] 24.63 24.64 3.25 2.79 19.40 29.66 1.00 1849 1614
## y1[8] 26.60 26.60 3.37 3.10 21.27 31.82 1.00 1827 1740
## y1[9] 28.59 28.60 3.71 3.21 22.69 34.35 1.00 1249 1409
## y1[10] 30.41 30.40 3.89 3.47 24.26 36.72 1.00 1075 1441
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -135.249
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 28 -18.8079 0.000961409 1.54484e-05 1 1 37
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -18.81
## b0 8.88
## b1 3.06
## s 3.18
## m0[1] 24.37
## m0[2] 20.62
## m0[3] 21.22
## m0[4] 31.87
## m0[5] 16.72
## m0[6] 30.70
## m0[7] 16.79
## m0[8] 17.92
## m0[9] 25.70
## y0[1] 25.42
## y0[2] 21.83
## y0[3] 20.08
## y0[4] 33.79
## y0[5] 17.84
## y0[6] 30.79
## y0[7] 13.15
## y0[8] 18.34
## y0[9] 23.62
## m1[1] 9.05
## m1[2] 11.98
## m1[3] 14.91
## m1[4] 17.85
## m1[5] 20.78
## m1[6] 23.71
## m1[7] 26.64
## m1[8] 29.57
## m1[9] 32.50
## m1[10] 35.43
## y1[1] 8.34
## y1[2] 8.50
## y1[3] 17.36
## y1[4] 20.74
## y1[5] 16.89
## y1[6] 25.68
## y1[7] 25.23
## y1[8] 29.82
## y1[9] 36.11
## y1[10] 37.43mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -19.50 -19.12 1.50 1.26 -22.58 -17.82 1.01 392 416
## b0 6.96 7.52 3.79 3.23 0.72 11.74 1.02 237 258
## b1 3.50 3.37 0.78 0.62 2.54 4.86 1.01 280 353
## s 4.49 4.11 1.58 1.21 2.69 7.56 1.01 566 682
## m0[1] 24.68 24.56 1.41 1.25 22.61 27.22 1.00 1137 1045
## m0[2] 20.40 20.46 1.45 1.25 18.05 22.64 1.00 427 637
## m0[3] 21.08 21.14 1.40 1.24 18.84 23.29 1.00 492 695
## m0[4] 33.27 32.79 2.68 2.28 29.86 38.48 1.01 575 484
## m0[5] 15.93 16.13 2.05 1.75 12.51 18.77 1.01 268 360
## m0[6] 31.93 31.50 2.43 2.10 28.81 36.61 1.01 640 508
## m0[7] 16.01 16.20 2.03 1.74 12.60 18.84 1.01 269 360
## m0[8] 17.30 17.44 1.83 1.56 14.24 19.95 1.01 290 385
## m0[9] 26.21 26.04 1.54 1.35 24.03 29.04 1.00 1207 871
## y0[1] 24.70 24.80 5.04 4.28 16.59 32.66 1.00 2019 1682
## y0[2] 20.33 20.46 5.07 4.49 12.08 28.26 1.00 1811 1642
## y0[3] 21.00 21.05 4.90 4.25 12.87 28.89 1.00 2100 1716
## y0[4] 33.31 32.92 5.34 4.80 25.42 42.14 1.00 1340 1389
## y0[5] 15.78 16.02 5.25 4.70 6.96 23.68 1.00 1249 1659
## y0[6] 31.96 31.49 5.46 4.80 24.08 41.18 1.00 1438 1354
## y0[7] 15.95 16.16 5.26 4.36 7.55 24.04 1.00 1024 1131
## y0[8] 17.20 17.40 4.92 4.53 8.81 24.70 1.00 1049 1289
## y0[9] 26.23 26.05 4.96 4.42 18.78 34.61 1.00 1789 1600
## m1[1] 7.16 7.72 3.75 3.20 0.98 11.89 1.02 237 252
## m1[2] 10.51 10.92 3.06 2.64 5.52 14.52 1.01 239 271
## m1[3] 13.87 14.14 2.41 2.06 9.85 17.13 1.01 249 284
## m1[4] 17.22 17.37 1.84 1.57 14.13 19.87 1.01 288 384
## m1[5] 20.58 20.65 1.44 1.24 18.28 22.82 1.00 442 644
## m1[6] 23.93 23.84 1.37 1.21 21.89 26.38 1.00 1014 1034
## m1[7] 27.28 27.05 1.67 1.48 24.95 30.47 1.01 1120 786
## m1[8] 30.64 30.23 2.20 1.90 27.79 34.89 1.01 729 558
## m1[9] 33.99 33.50 2.82 2.37 30.40 39.44 1.01 546 497
## m1[10] 37.35 36.75 3.50 2.89 32.90 44.06 1.01 459 447
## y1[1] 7.20 7.69 5.91 5.29 -2.96 15.64 1.00 550 841
## y1[2] 10.48 11.02 5.58 4.82 1.15 18.57 1.00 711 951
## y1[3] 13.99 14.27 5.42 4.77 5.01 21.89 1.00 859 1073
## y1[4] 17.22 17.51 5.03 4.44 9.02 24.86 1.00 1498 1452
## y1[5] 20.73 20.65 4.85 4.52 12.88 28.49 1.00 2050 1752
## y1[6] 23.83 23.80 4.80 4.13 16.50 31.64 1.00 1760 1554
## y1[7] 27.10 26.97 4.93 4.33 19.82 35.34 1.00 1638 1523
## y1[8] 30.68 30.37 5.27 4.48 22.47 39.55 1.00 1452 1151
## y1[9] 33.90 33.40 5.47 4.72 26.16 43.74 1.00 1328 891
## y1[10] 37.38 36.79 5.98 5.28 28.51 47.85 1.00 866 815
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
sensitivity/specificity analysis
ex14.stan
estimating sens and spec
data {
int N;
array[N] int x;
array[N] int y;
}
parameters {
real<lower=0,upper=1> p;
real<lower=0,upper=1> se;
real<lower=0,upper=1> sp;
}
model {
p~uniform(0,1);
se~uniform(0,1);
sp~uniform(0,1);
for (i in 1:N) {
y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
}
}
generated quantities {
real ppv;
real npv;
ppv=se*p/((se*p)+((1-p)*(1-sp)));
npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}n=20
data=list(N=n,
x=sample(0:1,n,replace=T),
y=sample(0:1,n,replace=T))
mdl=cmdstan_model('./ex14.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -26.386
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -11.803 0.0014423 8.29677e-05 1 1 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -11.80
## p 0.81
## se 0.23
## sp 0.43
## ppv 0.64
## npv 0.11mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.90 -17.59 1.29 1.10 -20.38 -16.44 1.00 731 895
## p 0.51 0.52 0.29 0.38 0.06 0.95 1.00 684 598
## se 0.27 0.25 0.11 0.12 0.10 0.47 1.00 2710 1570
## sp 0.44 0.44 0.15 0.16 0.19 0.70 1.00 2426 1502
## ppv 0.39 0.32 0.29 0.34 0.02 0.91 1.00 722 548
## npv 0.40 0.35 0.29 0.34 0.02 0.91 1.00 720 585
2x2 cross table
Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p01)
n11~B(n1.,p11)ex16-1.stan
data {
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p01;
real<lower=0,upper=1> p11;
}
model {
n01~binomial(n0,p01);
n11~binomial(n1,p11);
}
generated quantities {
real RR;
RR=p11/p01;
real OR;
OR=(p11/(1-p11))/(p01/(1-p01));
}n0=20
n01=rbinom(1,n0,0.3)
n1=20
n11=rbinom(1,n1,0.6)
data=list(n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -36.372
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 -22.217 0.00761187 0.000346553 1 1 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -22.22
## p01 0.15
## p11 0.55
## RR 3.67
## OR 6.92mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -26.62 -26.28 1.08 0.70 -28.81 -25.64 1.01 744 998
## p01 0.18 0.17 0.08 0.08 0.07 0.34 1.00 1182 823
## p11 0.54 0.55 0.10 0.10 0.37 0.71 1.01 1856 1223
## RR 3.83 3.17 2.89 1.56 1.45 8.36 1.01 1198 975
## OR 8.03 5.97 8.60 4.02 1.79 20.13 1.01 1315 1157
bimodal distribution, mixed normal distribution
n=100
y0=rnorm(n,0,1)
y1=rnorm(n,-5,1)
y2=rnorm(n,5,1)
y=sample(c(y0,y1,y2),n)
density(y) |> plot()
EM algorithm
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 3 components:
##
## log-likelihood n df BIC ICL
## -247 100 6 -521 -524
##
## Clustering table:
## 1 2 3
## 35 31 34rst$parameters## $pro
## [1] 0.346 0.308 0.346
##
## $mean
## 1 2 3
## -5.104 -0.269 5.189
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 3
##
## $variance$sigmasq
## [1] 0.953
##
##
## $Vinv
## NULLplot(rst)
ex17-1.stan
data {
int N;
vector[N] y;;
}
parameters {
simplex[3] h; //ratio of mix
real m1;
real m2;
real m3;
real<lower=0> s1;
real<lower=0> s2;
real<lower=0> s3;
}
model {
s1~cauchy(0,5);
s2~cauchy(0,5);
s3~cauchy(0,5);
for (i in 1:N) {
vector[3] p;
p[1]=log(h[1]) + normal_lpdf(y[i] | m1, s1);
p[2]=log(h[2]) + normal_lpdf(y[i] | m2, s2);
p[3]=log(h[3]) + normal_lpdf(y[i] | m3, s3);
target+=log_sum_exp(p);
}
}mdl=cmdstan_model('./ex17-1.stan')
data=list(N=n,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -707.624
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 39 -246.361 0.000153551 0.00108076 1 1 56
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -246.36
## h[1] 0.34
## h[2] 0.34
## h[3] 0.32
## m1 5.23
## m2 -5.13
## m3 -0.26
## s1 0.92
## s2 0.87
## s3 1.13mcmc=goMCMC(mdl,data,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 2100 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 2 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 3 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 4 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 2 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 3 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 2 finished in 1.7 seconds.
## Chain 3 finished in 1.7 seconds.
## Chain 4 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 4 finished in 1.7 seconds.
## Chain 1 Iteration: 101 / 2100 [ 4%] (Sampling)
## Chain 1 Iteration: 1100 / 2100 [ 52%] (Sampling)
## Chain 1 Iteration: 2100 / 2100 [100%] (Sampling)
## Chain 1 finished in 7.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 3.0 seconds.
## Total execution time: 7.1 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -255.32 -253.80 6.00 2.17 -275.94 -251.20 1.05 56 11
## h[1] 0.36 0.34 0.09 0.05 0.27 0.58 1.05 48 11
## h[2] 0.33 0.33 0.05 0.05 0.25 0.41 1.03 111 2246
## h[3] 0.31 0.33 0.09 0.05 0.04 0.41 1.09 29 11
## m1 2.47 5.12 4.45 0.35 -5.27 5.49 1.57 6 27
## m2 -0.13 -0.29 3.67 4.16 -5.29 5.36 2.41 4 25
## m3 -18.26 -4.92 78.70 1.08 -17.04 0.01 1.81 8 11
## s1 1.10 0.97 0.55 0.16 0.76 2.66 1.08 31 11
## s2 1.11 1.05 0.30 0.23 0.77 1.64 1.25 11 115
## s3 2.01 1.02 12.54 0.24 0.76 2.02 1.29 10 17
